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Abstract 

In this paper we propose a new mathematical model describing the effect of phosphoci- 
trate (PC) on sodium sulphate crystallization inside bricks. This model describes salt and 
water transport, and crystal formation in a one dimensional symmetry. This is the first 
study that takes into account mathematically the effects of inhibitors inside a porous stone. 

To this aim, we introduce two model parameters: the crystallization rate, which depends 
on the nucleation rate, and the specific volume of precipitated salt. These two parameters 
are determined by numerical calibration of our system model for both the treated and non 
treated case. 

Index terms — mathematical modelling, porous media, salt crystals, crystallization in¬ 
hibitors 


1 Introduction 

It is well known that one of the major causes of building degradation is the crystallization of 
salts into the porous matrix mmmm- Salt is present inside building stones as free ions: 
it can be a natural element of the material, created by reaction with atmospheric pollutants 
or introduced by water solutions penetrating into the porous matrix by capillarity [5j. The 
latter is the main mechanism leading to buildig damage and has received much attention from 
the scientific investigation [a 0 E, El M but remains not yet fully understood. Salt decay 
requires the simultaneous presence of soluble salts and water in the porous material, as well 
as appropriate environmental conditions. Indeed, it originates from salt-ions (e.g. chloride, 
nitrate, sulphate) that migrate while dissolved in liquid water which flows in the pore network of 
building materials. Liquid water may penetrate these materials by different processes, including 
hygroscopic moisture, penetration of rainwater (through, e.g. construction joints, damaged 
roofs, and cracks), dew point condensation, and rising damp. The latter is probably the most 
frequent and perhaps one of the most difficult sources of water to remove, when dealing with old 
buildings. Consider an initially dry porous stone (such as a masonry brick) that is wetted by a 
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salt water solution. During the wetting phase, water fills up the stone bringing the dissolved salt 
present in the outside environment. If the stone is in contact with ambient air, water molecules 
are exchanged with the environment by evaporation thus starting a drying phase; the rate of 
dehydration depends on the relative humidity of the atmosphere. At this point, salt content in 
water increases and solution may become supersaturated. Once a high degree of supersaturation 
is reached, salt starts crystallizing: if crystals are formed inside the porous matrix we talk of 
subflorescence or cryptoflorescence ; if crystallization takes place on the exterior boundaries of 
the stone we talk of efflorescence. Subflorescence causes the formation of large crystals into 
the pores: once the presle exerted by these crystals exceeds the tensile strength of the porous 
matrix, it can lead to widespread loss of surface, e.g. exfoliation, detachments. The occurrence of 
efflorescence or subflorescence (cryptoflorescence) depends on several factors including salt type 
and concentration, microclimate, evaporation rate |6], substrate porosity characteristics mma 
011 and surface tension and viscosity of the solution [El ESI EH- The in-pore crystallization 
causes a reduction of the pore volume, breaking the liquid network and delaying water transport. 
Since pore clogging affects the location and quantity of crystals, it might have implications 
for stress development and deterioration of the material fl8 j. Common constructions contain 
different kind of salts such as chlorides, sulphates, nitrates, and carbonates, with their own 
solubility, crystalline structure and crystallization properties. Among these, sodium sulphate 
is probably one of the most complex and damaging salt types involved in salt decay processes. 
Indeed, this salt has three different phases of crystallization at various microclimate conditions, 
can easily supersaturate and has a solubility which is highly temperature dependent HSUU]. 
Both crystallization and hydration transformations in sodium sulphate, resulting in significant 
volumetric changes, have been blamed for the destructive mode of action of this salt {21]. 

One way to prevent the stone breakage is to treat the porous material with a substance 
that inhibit subflorescence: these crystallization inhibitors reduce the pressure associated with 
the growing crystals trying to keep it below the breakage modulus of the substrate. The or¬ 
ganic as well as inorganic ion and molecule additives alter the surface properties of the crystals 
which lead to changes in nucleation, growth, and thereby changes in the shape of the crystals 
as well as in their agglomeration/dispersion behaviour. Examples of well-known additives with 
extended technological and industrial uses are the families of (poly)phosphates, carboxylates, 
polyacrylic acid derivatives, and benzotriazoles [221122 23 ] • These additives are widely used 
as scale-inhibitors to prevent undesired effects associated with sparingly soluble salts (e.g. sul¬ 
phates, carbonates) precipitating in oil extraction pipelines [25] industrial boilers, heat exchang¬ 
ers, house appliances or water pipes [261 127] and others. The effectiveness of a given inhibitor 
depends on many variables: salt type, pore structure properties of the substrate, application 
methodology, the composition of the inhibiting solution to cite a few. Hence, a given modifier 
has to be evaluated for each stone and for each salt [28]. On the other hand, adding a crystal 
inhibitor does not affect surface tension nor contact angle of the wetting liquid, since there have 
not been observed any significant effect on solution transport |29| . Although the effectiveness 
of some salt crystallization inhibitors in bulk solution has been proved, the possibility of using 
these products for the prevention of salt decay in building materials is still controversial because 
it is not clear how these inhibitors act. However, experiments suggest two possible mechanisms 
[301E21 ED M- nucleation delay enhances salt transport toward the surface, thus increasing 
efflorescence; crystal habit modification by absorption on specific faces of a growing crystal that 
decreases crystal growth rate. Another matter of discussion is the fact that crystal reduction 
would result in higher supersaturated solutions. It has been speculated, but not actually ob¬ 
served, that in this case the inhibitor may promote salt precipitation at higher supersaturation 
levels and, hence, the quick formation of large crystals. Therefore, a modifier would eventually 
increase the crystal pressure and the risk of damage instead of reducing it [33|. Our group has 
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undertaken a broad research project |34l 135 , 36] focusing on the effects of environment-friendly, 
non-invasive inhibitor systems on saline solutions percolating and crystallizing in a porous media 
following evaporation, in order to develop a sound methodology suitable for addressing the con¬ 
servation needs of different salt-weathered sites. Our attention has been particularly focused on 
the crystallization inhibition properties of functionalized polycarboxylates (i.e. maleate, citrate, 
phosphocitrate, tartrate), with an emphasis on the phosphorylated family members. Indeed we 
have demonstrated that phosphocitrate (PC) has been revealed to be one of the most promising 
inhibitors, because of its effectiveness in controlling the crystallization of different salts (i.e. 
sodium sulphate, sodium chloride, sodium nitrate, calcium carbonate) and salt mixtures in a 
wide range of porous materials and in various ambient conditions. 

In this work we developed a mathematical model describing the effect of phosphocitrate 
(PC) on sodium sulphate crystallization inside a brick’s porous matrix. There are plenty of 
mathematical models describing salt crystallization in porous stone. They consists of 3D mul¬ 
tiphase systems of equations for heat and mass transport with various degree of complexity. 
Some models might also couple the governing equations with other effects: osmosis, stress ten¬ 
sor deformations and latent heat release due to salt crystals formation [33 EM SDI. For the 
present study we have developed a simple mathematical model of salt and water transport and 
crystal formation. In fact, we limit our research to the considerations on few available data, 
which can be obtained using simple laboratory equipments, and so it would not have made 
sense to include further effects. Moreover, since the experiments were carried out in laboratory 
at constant temperature, we did not consider directly temperature variations; we just included 
evaporation rate into the porous stone simply by defining an appropriate sink term in the water 
balance equation. Actually, this work is a preliminary study to describe mathematically the 
effects of inhibitors inside a porous stone: to our knowledge, this is the first attempt to develop 
a mathematical model for the effects of crystallization modifiers. As we shall see, we identified 
two model parameters that will be crucial for the appropriate description of an inhibitor: 

K s the crystallization rate taking into account the nucleation rate; 

7 the specific volume of precipitated salt, describing the crystal habit modification. 

These two parameters will be determined by the numerical calibration of our model - i.e. by 
comparing our numerical results with the available experimental data - for both the treated 
and non treated case. The remain of our paper is organized as follows: the second section will 
describe the materials considered and the experiments performed; the third and fourth sections 
will introduce the mathematical model and describe the numerical scheme applied to solve the 
system equations; in section five we will describe our results. The paper ends up with few 
conclusions. 

2 Materials and Methods 

In this section, we will introduce the experimental settings we are going to consider m- Com¬ 
mercially produced brick is tested. Bulk density p v was determined by weighing and measuring 
of dimensions of dry prismatic samples. The matrix density p ma t was measured by helium 
pycnometer. The porosity no [%] was calculated according to the equation 

n 0 — 100 • (1 p v /Pmat) • ( 1 ) 

The porosity determined in this way is 28.51% ± 0.04%. 

Pore size distribution was determined by mercury intrusion porosimetry (MIP) by Carlo 
Erba instrument on a about lg of material. All experiments were performed in air conditioned 
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laboratory at 25±2°C and 30±5% RH. Table[l]shows the pore size distribution of the considered 
brick. 


Pore Radius 
Interval (/irn) 

0 .001-0.01 

0 .01-0.1 

0.1-1 

1-10 

Distribution (%) 

2.8 

9.3 

42.2 

45.7 


Table 1: Pore size distribution in the brick under consideration. 


2.1 Experiment 1: brick’s capillary absorption and drying test in pure water. 

This set of experiments were conducted, according to standard UNI EN 1925 (Determination of 
water absorption coefficient by capillarity) and NORMAL 29/88 (Drying Behaviour), without 
the presence of salt. It will serve as a control sample to test transport properties of the materials 
under study. The brick specimen has the form of a cube of side 5 cm, is positioned in a bucket 
containing water and immersed for 3 mm in height. The water absorption for capillarity, 
expressed in g/cm 2 , is defined as the quantity of water absorbed by the specimen having the 
base surface in contact with water as a function of time t, with room temperature and pressure. 
At different time intervals the specimen is taken and tamponed only on the wet surface and 
then weighted until the variation in the quantity of absorbed water between two consecutive 
measurements, for a 24 hours interval, is less than 1% of the water mass. The determination 
of the quantity of water absorbed by the specimen per time unit is given by W = 
expressed in g/cm 2 , where W is the quantity of water absorbed (expressed in g) and S = 25 
cm 2 is the surface of the specimen in contact with the porous frame. The experiment is applied 
to a number of specimen and then the average of the time dependent values W obtained for the 
different specimens is computed. Finally we get the averaged quantity with tf. the time 

instants expressed in s 1 ^ 2 . 

2.2 Experiment 2: brick’s capillary absorption and drying test in a salt 
saturated water solution. 

Both in untreated and treated brick’s samples with PC the water and salt concentration profiles 
were determined experimentally using prismatic specimen 2 x 2.5 x 12 cm positioned vertically 
in a bucket containing a salt water solution of Na, 2 SO± (99.5 g/L) (see Fig. [3]). In order to 
determine the concentration profiles the specimens were cut into 4 pieces with similar dimension 
and re-assembled sealing the lateral sides with epoxy resin; in this way only the top side of the 
brick is in contact with ambient air. On the other hand, the immersed part of the specimen 
is pervious and liquid can flow through the lateral side. The insulated specimens were dried 
at 65 ± 2°C to the constant mass. When the solution in the bucket is totally absorbed by the 
specimens, the water content was obtained as difference of the mass of the saturated specimens 
and of the sample’s mass after drying at 110 ± 2°C to the constant weight. The concentration 
of sulphates in the dried samples was determined as follows: the samples were placed in plastic 
container, 200 mL of boiling water was added and the container was sealed. This procedure 
was repeated every day for 1 week. Then the dry samples were weighted and the concentration 
of sulphates was calcutated. 
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3 The mathematical model 


Here we want to introduce a model of coupled water and sulphate transport taking into account 
not only the influence of water flow on salt transport but also the effect of bound sulphates on 
pore walls, and the effects of porosity changes (due to the salt bonding) on moisture transport. 
Regarding the mathematical domain, a reasonable assumption is to consider a one dimensional 
geometry since the domain is sealed on its lateral side; hence, flow is predominantly vertical. 
We denote by n the porosity, i.e. the fraction of volume occupied by voids, and we denote 
the fraction of volume occupied by the liquid and by the gas (composing the fluid) within the 
representative element of volume, respectively by 81 and 6 g . The following relation holds: 

n = 61 + 9 g. (2) 

The mass balance equation for a liquid of density pi reads as: 

=/(^) ( 3 ) 

where q is the water flux into the porous matrix and f ( 9 1) is the evaporation rate inside the 

specimen. Both q and f(6j) will be specified later on. 

Let us denote by <7 the concentration of free ions in water and with c s the density of bound 
salt, the mass balance equation for salt dissolved in water is given by: 

c)c 

d t (0ia) + d z (ci<\) = Dd z (0id z Ci ) - ( 4 ) 

at 

where D is the salt diffusion coefficient, while the sink term on the right hand side takes into 
account the crystal formation into the porous matrix. In this work, we assume that crystal 
growth depends on the following properties: the concentration of salt dissolved in the liquid, 
the fraction 9 g and the degree of supersaturation. If we indicate the supersaturation level with 
c we have: 

f)r — 

= K s c i e 2 g +K(c i -c) + e l . (5) 

with K s and K two crystallization coefficient and (•)+ is the positive part function (or the second 
term is active only when salt saturation into the liquid exceeds the supersaturation level). The 
term 9 g on the right hand side is raised to power two in order to capture the following fact: the 
higher the water content, the smaller the crystallization into the pores. The power two simply 
slows down the crystal formation in saturated regions. The second term on the right hand side 
has been defined for the sake of completeness; in fact, in our experiments and in the subsequent 
simulations, salt supersaturation has never been exceeded and term K has not been determined. 

Since the overall porostity changes as the salts growth into the porous material, the following 
equation holds: 

n(t) = n 0 - jc s , ( 6 ) 

with 7 the specific volume of sulphate crystal. 

3.1 Darcy’s law 

Water flow into a porous medium is given by the well known Darcy’s law [42, 33] : 

q = ) {d z P c (s) - pig) (7) 

m W / 
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with P c = P c (9i/n ) the capillary pressure, k the permeability of the porous matrix, hi the 
viscosity of the fluid, the term (n/no) 2 is a shape factor for the influence of the porosity variation 
to the water flux and s = Qi/n. 

Capillary pressure is usually given as a function of water saturation and is defined through 
a state equation. In literature, one can find capillary pressure state functions for several ap¬ 
plications; in building materials, however, despite the number of experimental study, there is 
not a relation correlating capillary pressure with moisture content into the porous matrix. To 
overcome this problem, we will approximate Darcy’s law through a polynomial function with 
some free parameters that will be found through model calibration. Thus, we proceed as in |44j . 
First of all, since the dimensions of the brick are small, gravity effects can be safely disregarded 
from ([7]). Then we introduce function B such that 

d z B = J±ld z P c {-). 

Pi 


We know that P c (s ) is a decreasing function of s = 81 /n < 1 and vanishes whenever the medium 
is completely saturated, i.e. 9i = n. On the other hand, permeability k = k(s) is a non-negative 
increasing function of s and it is bounded from above by its value at saturation. Taking into 
account these observations, the first derivative of function B with respect to s = Oi/n can be 
given by the ansatz 


B (s) = max 


4c 


(a 


s)(s - 1), 0 


( 8 ) 


with a such that k(a ) = 0. Constants a and c are physical properties of the porous material 
involved and will be determined later on. The quantity a ■ n is the minimum value for saturation 
ensuring the hydraulic continuity - i.e. water transport through the porous medium. On the 
other hand c has the dimensions of a diffusivity. The term 4c/(l — a) 2 is chosen so that 
maxl-B^s)} = c. Integrating B'(s ) we obtain the following expression (see Fig. [I]): 


! c ((r=f) (3a - 1 - 2s) + (1 - a)l, if s G [a, 1], 

B(s) = < n V rn n J 

K 0, if s G [0, a), 

B{ 1) = |c(l — a), if s > 1. 

Summing up, Darcy’s law can be expressed as follows: 


q = 




3.2 Water evaporation 

Once water content decreases below the quantity an, the hydraulic continuity is broken and 
fluid trasport is no longer ensured. Since drying experiments end up with a completely dry 
stone, we added a sink term in the water balance equation © to take into account the effect 
of evaporation inside the porous matrix. In our mathematical model, we made the simplifying 
assumption that evaporation is maximum when moisture content is below the value an and 
decreases quickly as the porous medium becomes saturated: thus liquid flow and evaporation 
acts at almost separated stages (one is strong while the other is weak and viceversa). This is 
reasonable since in our controlled experimental setting temperature is constant and does not 
play a significant role. 

We defined the evaporation rate as follows 

f(0i) = — PiKt9iH £ (9i) (9) 
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x 10' 




Figure 1: Graph of the functions B' a c {9) (above) and B a ^ c (9) (below) for s E [a, 1] for the choice 
a = 0.219 and c = 9.87- 10~ 4 . 


with Kt a (temperature dependent) constant and 7~L £ is defined as follows: 


H s (9) 


1 

an+e -1 n , (an+e)(l-an) 
£ £ 
an+e n _ n(e+an) 
n(a— l)+e n(a— l)+e 


if 0 < 6 < an, 
if an < 9 < an + e 
if x > an + e. 


see Fig. [2} In our simulations, we took e = 0.25 an. 


( 10 ) 



Figure 2: Profile of the function 'H £ (9), with e = 0.25 an. 
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Description 

Units 

Value 

Ref. 

hi 

Brick’s height in the experiment 1 

cm 

4.7 



h -2 

Brick’s height immersed in the solution 

cm 

0.3 



hi 

Brick’s height in the experiment 2 

cm 

11.7 



n 0 

Porosity of the unperturbed material 

- 

0.2851 

Eq. 

2 


D 

Diffusivity of Na 2 SC >4 

cm 2 /s 

1.230 x 10 “ 5 

|4ol sect. 

6 . 2 ] 

Pi 

Density of water 

g/cm 3 

1 

m_ 

9i 

Moisture content of the ambient air 

g/cm 3 

6.254 x 10 ~ 2 

Eq. 

27 


c 

Saturated concentration in water of sodium sulphate 

g/cm 3 

0.4399 

Eq. 

5 


Ci 

Concentration in water of sodium sulphate 

g/cm 3 

9.95 x 10 -2 

Eq. 

18 



Table 2: Parameters of the problem. 



Description 

Units 

a 

physical property of the porous matrix 

- 

c 

physical property of the porous matrix 

cm 2 /s 

7 

Specific volume of crystal 

cm 3 /g 

Ki 

Exchange coefficient 

cm/s 

K s 

Crystallization rate coefficient 

s _1 

K 

Growth rate of hydrated crystals 

s -1 

Aj 

Evaporation rate 

s -1 

a 

Evaporation exponent 

- 


Table 3: Model coefficient to be calibrated 

3.3 The Complete Mathematical model 

Summing up, the mathematical model we are going to consider is the following: 

d t 9i = d z ((£) 2 &£( 0 iAO) - Kt9,H £ (9,X 

< d t (9ici ) = d z (ci d z B(0i/n) + 6iDd z c^j - (11) 

n{t) = no — 7 c s , _ 

„ °gt = K s Ci(n - 9i ) 2 +I\{ci - c)+6i. 

Table [ 2 ] shows the known parameters of the problem. Since some coefficients are unknown, 
we will calibrate the model versus experimental data. The obtained values will give some insight 
about the action of the inhibitor in the crystallization process. Table [3] lists the coefficients to 
be determined. 

3.4 Boundary Conditions 

For each experiment we will describe the initial and boundary conditions to apply to model 
In some cases, we are even able to simplify the model equations. 
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3.4.1 Experiment 1: pure water 


The immersed part of the brick (for — h 2 < z < 0) is pervious to later water flow and we assume 
that it is initially saturated. In this way we can simply confine ourselves to mathematically 
describe the domain 0 < z < h±. Moreover, since there is no salt, our mathematical model 
reduces considerably; indeed, we can only retain the water continuity equation ([3]), that in this 
setting, is given by: 

d t 6 i = d zz B - K T 6 {H e { 6 l ). (12) 


Given the absence of salt, porosity will remain constant and, thus, will not affect water flow. 
Equation (12) has to be coupled with reasonable initial and boundary conditions. For the 
experiment of imbibition, we assume the conditions 


r <9i(z, 0) = 0, 

\ 6 i( 0 , t) = n 0 , 


(13) 


that is, the sample is initially dry while its botton side is always saturated. To reproduce the 
loss of water at the upper boundary z = h\ due to evaporation, we derive 6 i(hi,t) from the 
following relations: 

r d z B = km - Oil*- 1 - di), if o t > e h 

\ 61 = 0 i, otherwise. 

In the above conditions, 9[ is the moisture content of the ambient air (assumed constant) while 
K[ is the exchange coefficient with the environment. The exponent a > 1 takes into account that 
water evaporation from the top of the doamin depends non-linearly on the difference between 
the quantity of water within the specimen and the value 6 [. 

Once the imbibition stage is terminated, we stop the simulation and switch to another 
settings to deal with drying. In this case we consider the whole domain [— h 2 ,h\\, since we 
do not add water at the bottom of the specimen. The other changes regard the initial and 
boundary conditions. If we denote by t s the final time of imbibition and with 0fi n (z) = 9i(z, t s ) 
the value of Q\ after imbibition, the initial condition for the new setting is given by 


( 6i(z, 0) = 0 fi n (z), for z e [0, hi], 
1 0 i(z, 0) = n 0 , for z e [~h 2 , 0 \ 


meaning that the initial water content is the final value obtained for the imbibition test. More¬ 
over, at z = — h 2 we impose a no-flux boundary condition: 


d z 9i(-h 2 ,t) = 0, 


while at z = hi we retain condition (14) again. 


(16) 


3.4.2 Experiment 2: salt saturated water solution 

Experiments with salt solution were performed on bricks with height 12 cm. As above, we will 
consider during imbibition that the first three millimiters are submerged with water, thus we 
confine ourselves to the domain [0, /13], while during evaporation, to the domain [—h 2 ,h 3 ] (see 
table [ 2 ]) . 

For imbibition, we assume the initial conditions for the system 

< 


c s (z, 0) = 0, 
Ci(z, 0) = q, 
0t(z, o) = e h 
n(z, 0) = no- 


(17) 
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Figure 3: Setup of experiment 2 as described in Section 2.2 


As boundary conditions for t G [0,t s ], we impose for the ion content, at z = 0, the salt 
concentration q of the solution used in the experiment: 

Cj(0, t) = Ci (18) 

with Ci the actual concentration of sodium sulphate in water and a saturation condition for 
the water content 


= n(0,t). 

At the top boundary z = h 3 , we impose: 


(19) 


d z Ci{h 3 ,t ) = 0, 


i.e. zero ion flux through the upper brick boundary and condition (14). 
For the drying phase, we assume the initial conditions: 


c s (z,0) = c s (z,t s ), 
Ci(z, 0) = Ci(z,t s ), 
Qi{z, 0 ) = 0i(z,t s ), 

„ n(z, 0 ) = n(z,t s ), 


( 20 ) 


( 21 ) 


with z G [0, h 3 ] and for the immersed part, corresponding to z G [—h-2,0], of the specimen we 
set: 


c s (z,0) = 0 , 

Ci(z, 0 ) = Ci, 

0i(z, 0 ) = n 0 , 
„ n(z, 0 ) = n 0 . 


( 22 ) 


From now on we consider separately the four bricks composing the specimen. To this aim we 
define as h\ the height of the broken brick and the points 


hi = with <n = 0 for i = 1, 2, 3,4. 
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The *-th brick is then parametrized as the interval [a,i,bi\ for i = 1,2, 3,4. Then as boundary 
conditions we impose at the bottom z = Oj, zero ion flux through the lower brick boundary 

d z Ci(cii,t) = 0, * = 1,2,3,4 (23) 

and as a boundary condition for 9\ reproducing the loss of water at the lower boundary we 
assume 

e l (a i ,t) = e u i= 1,2, 3,4. (24) 

At the upper boundary z = 6* we assume the conditions 

cU(M) = 0, * = 1,2, 3,4 (25) 


and 

and we put 9i 
the oven. 


e l {b u t) = e u i = 1,2,3 ,4 


(26) 


0 in both conditions (24) and (26), in order to reproduce the situation inside 


3.4.3 Calculation of parameter 6 [ 


Since we do not have any measurements of the relative humidity of the ambient air surrounding 
the sample, we set the value of the moisture content in the environment using the value of 
the average quantity of water within the brick measured in the imbibition-drying experiment 
with the sole water. In particular, using the measured average value Q s (quantity of water at 
saturation of the specimen) and the average value Q d = Qf m — Q s (loss of water at the end of 
the drying the experiment) we compute the final quantity of water Qf in = Q s + Q d = 0.31274 
g/cm'? and then we get: 


Si 



0.06254. 


(27) 


4 Numerical approximation 


Here we propose a numerical scheme for the model ©• We mesh the interval [0, h\ with a step 
A z = jj and we denote 

At 

A = — ,Zj =jAz,j = 1 

We also set w k = w(zj,tk ) the approximation of the function w at the height Zj and at the time 
tk- As showed in m The simplest and consistent approximation of d z (r(z)d z w ) by means of 
Taylor expansions is the following first order approximation: 


Aj(r,w) : = 


0 rj + r j+ i)(w j+ i - wj) - (rj-i + rj){wj - wj- 1 ) 
2A^ 2 


(28) 


From now on, we will omit for simplicity the subscript l of 6 . Then, the discretization in explicit 
form the first equation of the model is: 

^fc+l _ gk 

j At j = Aj((n fc /n 0 ) 2 , B k ) - A tK T H £ (e k )6 k+1 , (29) 

Now, if we consider the velocity field computed in the equation pi)]) and we set it as V = 
(■ n/no) 2 d z B(9/n ), we can rewrite the second equation of the system (11) as: 


dt(0ci) - d z (ciV) = d z (D6d z Ci ) - K s a(n - 0) -K(a - c)+0. 


(30) 
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We can assume: 


47 =i — 


n 0 


2 B 


°J+1 


7+1 


- B 


°U 


"j-1 


2Az 


, for j = l,...,JV-l, 


(31) 


with the boundary values set as follows: 


T/fc _ 
1 /q — 


0 , for the imbibition phase, 


- ( ^ ) I<i{0i - 0 k ), for the drying phase, 


and 


n 


Va) = ( ) KAOi — 9j), for both phases. 

no 1 


(32) 


(33) 


Therefore, an explicit and monotonic scheme for (30) reads as: 


(0ci)* + L - ( 0*) 

At 


k t rk „k _ t rk „k 

' A —1— 1 o _L 1 ’ n _ 


v j+lH,j+l 


j -1 »J-1 


+ 


2 Az 

i^+ii^j+i - 2 \ v Mi + 


2A z 

Jz 


(34) 


+ A J (Zl^,cf)-A>f J (n J fc -0|) 

— K( c ij — c)+ 9j, 


which is convergent under the CFL condition 

At < 


inf 9jAz 2 


Dn o + sup\V\Az + (K s + K)noAz 2 


We observe that the CFL may become very restrictive during the drying phase, since 9j tends 
to zero. For this reason we simulated separately the two phases (imbibition and drying) using 
two different lower bounds for the CFL taking into account the evolution of 6 j in the two cases. 
Then, using the Euler’s method for approximation of the third equation in ( |11[ ), we can write 
the discretized problem as: 


0 k+l = Q) + At Aj{(n k fno) 2 , B k ) - A tK T n e (9 k )9 k+ \ 
j = 1 __ 

= c k s j + A t[K s c k j (n k - 9 k ) +K(<* d - c)+ 0% 
j = 0,...,N 

= n 0 - 7cjJ 1 , j = 0, • • •, N 

i_ / ok k A l^d< ?+1 -2|v/|<,+|L/AI< ? ^ 

: +! 1 7/ c i,3 ^ A 2 


J 

r k +1 

"‘ +I 

A+l 


(35) 


gfe+i 
3 


+AtAj(D0 k ,c k ) + \ 


yk k _yk r k 
k „k\ I \ v j+l c i,j + l v j—l c i,j—l 


x 3\ ± ^ w i '-'i ) ' 2 

At[K s cUn k -9 k )+K(c k ,-c) + 0 k 


3 3 - 

j = 1,...,AT- 1, 


l}' 


with suitable boundary conditions described in the next subsections. In particular, for the 
first equation of the scheme we have: 


9 k+1 = C(0j + At Aj((n fc /n 0 )^, B k )) 


(36) 
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with 


1 


C = 


1 +a tK T nMy 


Note that the scheme in the last equation of (35) may become degenerate if 9 k+1 is null, 


thus we put into the numerical algorithm a threshola for 9 in order to avoid this possibility. 

4.1 Boundary conditions for the imbibition phase 

At the bottom boundary of the brick, we assume the condition for the ion content according to 


the concentration value of the experiment (18), which reads as 


c^+i — r- 
L i, 0 — 


(37) 


and the condition (47). At the top boundary of the brick, we impose the zero ion flux condition 


(20) for the ion content, discretized with a second order approximation: 

Jc+i _ £ r k+i _ 1 fc+i 
c i,N - £ L i,N-l 3 C i,N- 2 ' 


(38) 


Let us now consider the discretization of the condition (14), reproducing the exchange with the 


environment. Note that in the case of the experiment 1 with sole water in the condition (14) 
we have to replace n !■ with the constant value no- 


At the node zjy we need to solve the equation 

3 


2Az 


B{9 k + 1 /n k + 1 ) + K x % - 9 k N +1 \ a -\9i - 9 k N +1 ) = 
^B(9 l p~_} l /n l p'_} 1 ) — B{9 k ^ 2 /n k j^\) 


(39) 


2 Az 


with the function to be inverted 


9 i W = - Ki\9 t - 9\ a -\9 t - 9). 


2Az 


The invertibility condition is 


9\ = 


d e B{9/n) + K l a\9i-9\ a ~ 1 >0 


2nAz 

on a compact set, with 

doB(s = 0/n) = 


(40) 


J / a+ 1 _6 _/ 9 \ J _ a_ 

I no «o \ n °) n ° n ° 

0 , elsewhere . 


if s E [a, 1], 


(41) 


Note that the condition dgB(6/n) > 0 is always satisfied for 9 E [o • n,n\, so that (40) holds. 

ton’s 

(42) 


Therefore, at the upper boundary of the brick we need to solve, using for example with Newton’s 
method: 


() k + l - n~ l 

6 n ~ 9 1 


4 J B(^+ 1 i/n^L 1 1 ) B(9 k ^ 2 /n k ^ 2 ) 


2Az 
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4.2 Boundary conditions for the drying phase 


In order to model the loss of water, we use the zero ion flux at the bottom of the brick, discretized 
with a second order approximation as 


k+i f fc+i _ 1 fc+i 
L i,o 3 °*’! 3 S;,2 > 


(43) 


and condition (38) at the top boundary. Let us now consider the discretization of the conditions 
(24) and (26), reproducing the situation of the specimen inside the oven we set at the lower 

8 k o +1 = 0 U (44) 


boundary: 


and analogously at the upper boundary: 


ofc+l _ n 


(45) 


with 61 = 0 . 


5 Numerical Results and comparison with experimental data 

5.1 Calibration of parameters a,c,Ki, Kt,&- 

Now we describe the calibration procedure to determine a, c, Ki, Kt and a for both the phases 
of imbibition and evaporation of water in the brick using the experimental data of experiment 
1 . 

We need to compute the total quantity of water absorbed and lost by the brick at time t k given 
by: 


rh i 


pid(z,t k )dz, 


(46) 


thus we need to solve problem (12). We compute 9(z,t k ) numerically with the forward-central 
approximation scheme 

At 


9 i +1 = d * + ^ 


(B a , c (6j +1 /n 0 ) - 2B atC (&j/no) + ^a,c(^-i/n 0 )) 


with the boundary condition at the top boundary (42) under the CFL condition 

At n 0 _ no 

Az 2 ~ 2 d z B a>c ~ 2c’ 


with 9 = 9(zj,t k ), zj = jAz,j = 0^ , {4}fc=i,...,JV meos - At the bottom boundary 


we use the imbibition condition 


^' +1 = n£ +1 


'o ~ '"0 ( 47 ) 

and Neumann condition 9~(0,t) = 0 of null flux, only for the drying phase, that numerically 
results to be 

A . . 1 

(48) 


afc+1 4: nk+1 ^ flk+1 

Vn — —V 1 — —Vo 


3 1 3 

Let us define t s the saturation time at the end of the imbibition phase and Q s the corresponding 
value. Then we compute the approximated values of the quantity of water in the brick Q% um 
at time t k as follows. With the trapezoidal rule we compute the integral (46): 


JV-l 


or™ = p— («o +2 £ 

3 = 1 
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Quantities 

Value 

Dimensions 

a 

0.21904 

- 

c 

9.8073 x 10 “ 4 

2 —1 
cm s 

Ki 

3 x 10 -5 

s -1 

Kt 

3.2 x 10 " 7 

s -1 

a 

0.9 

- 


Table 4: Results of the calibration for the imbibition and drying stages without salts. The 
overall error is about 7%. 


u 


* data 

• fitting values 


0 200 400 600 800 1000 1200 1400 1600 

time [a 1 '' 2 ] 


Figure 4: Data fitting result: comparison between data points and fitting values obtained for 
a = 0.21904, c = 9.8073 X 10" 4 , K t = 10" 5 , K r = 3.2 x 10" 7 , a = 0.9. 


in order to compare the numerical quantity of water to experimental data Qk at time tk- The 
error to be minimized is then defined as 


E(a,c,Ki,K T ,a) 


1 'yT \Ql um -Q k \ 
Nmeas \Q\ 


with Q the average value among data. The calibration procedure has been carried out in 
MATLAB© applying the simulated annealing method. The computational time for a single 
simulation with fixed parameters takes 900 seconds on an Intel(R) Core(TM) i7-3630QM CPU 
2.4 GHz. Table i lists the results obtained within an error of about 7%. 

Figure [4] shows the comparison between measured data and numerical simulations after 
calibration. 


5.1.1 Calibration of constants K s and 7 


As described in Section 2.2 for experiment 2 , the bricks were first broken in four pieces with 
similar dimensions, both for the treated and non treated cases; for any brick, we measured its 
salt content. 

In order to determine constant K s and 7 we need to define an appropriate functional to be 
minimized. We proceed as follows. First we define the average quantity of salt in i-th brick as: 

rbi 


A 

Bi 


c s (z, t)dz = 


1 

/if 


c s (z , t)dz for 7 = 1,2, 3,4. 


(49) 
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where t is a sufficiently long time when we can assume that the water is completely evaporated. 
Here Ai and Bi represent the cross section and the volume of brick i, respectively. 

If we denote by qV’ um the average quantity of salt in i-th brick obtained discretizing formula 


(49) with the trapezoidal quadrature rule, the values of K s and 7 can then be found solving the 
following minimization problem 


with a the average salt 


1 4 Lnum r , 1 

■ 1 m ~ Qi 

mm - > - — 

Ks,l 4 f—f |<?| 


(50) 


i =1 



0.15 


S °- 1 i. 

a 

| 0.05 - 


o'- 1 - 1 - 1 - 1 - L - 

0 2 4 6 8 10 

height of the specimen emeiging from tic solution cm] 


Figure 5: Experiment 2. Imbibition phase in the salty solution: profile of 61 and 7 depicted at 
the final time of the experiment T = 1128 h, with K s = 4.1 • 10 ~ 5 s -1 and 7 = 0.6. 


For the case without inhibitor, the profiles of the quantities obtained numerically at the end 
of the imbibition experiment (47 days) for the not treated bricks, indicated by NTi, i = 1,2, 3,4, 
are depicted in Figg. [5] and [6] As expected, the quantity of water in the brick 9i is a decreasing 
function of the height of the brick, since the top of the brick is interested by water exchange with 
the exterior. The graphs of the same quantitites at the end of the drying phase are depicted 
in Figg. [7] and [8j We observe that the amount of bound salts is, as expected, an increasing 
function of the height of the brick, since crystals mostly form where the quantity of water is 
lower. The calibration procedure gives the following result: we obtain an error of about 11.6% 
for the values K s = 4.1 • 10 -5 s -1 and 7 = 0.6 cm 3 g _1 . In Table [ 5 ] we report the comparison 
between measured data and numerical values obtained using the parameters deriving from the 
calibration procedure. 

For the experiment of the bricks treated with PC-10~ 6 M, at the end of the calibration proce¬ 
dure we obtain an error of about 13.7% for the values K s = 6 • 10~ 5 s -1 and 7 = 0.53 cm 3 g _1 . 
In Table [6] we report the comparison between measured data and numerical values obtained 
using the parameters deriving from the calibration procedure for the four bricks, indicated by 
PCi, i = 1,2,3,4. In Fig. [ 9 ] we depicted the profile of c s for the not treated case (NT) and in 
presence of PC-10 -6 M (PC). As observed experimentally, the amount of salt crystals is higher 
in the case of the treatment with the crystallization modifier. 


16 







0.03 




Figure 6 : Experiment 2. Imbibition phase in the salty solution: profile of c s and n depicted at 
the final time of the experiment T = 1128 h. with K s = 4.1 • 10 ~' 5 s -1 and 7 = 0.6. 
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Figure 7: Experiment 2. Drying phase: profile of 61 and Cj, with K s = 4.1-10 5 s 1 and 7 = 0.6. 
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quantity 

NTl 

NT2 

NT3 

NT4 

Qi 

14.62 

17.18 

17.74 

30.18 

n num 

12.21 

17.88 

22.52 

30.10 


Table 5: Salt content in any small brick in the not treated case (NT). We reported the measured 
salt content qi and the numerical values qV- um expressed in mg/cm 3 . 


o.i 
0.08 
f 0.06 
I 0.0'' 
* 0 . 02 ' 


- NT1 
NT2 
-NT3 
■ NT4 


^ £ 5 5 o o o ° ~ < 


0 0.5 1 1.5 2 

2.5 


—1— NTl 

- NT2- 

-«-NT3 

—•— NT4- 




0.5 


1 1.5 2 

height of the specimen [crn| 


2.5 


Figure 8 : Experiment 2. Drying phase: profile of c s and n, with K s = 4.1-10 5 s 1 and 7 = 0.6. 


The computational time for a single simulation with fixed parameters both for the treated 
and the not-treated case takes 2240 seconds on an Intel(R) Core(TM) i7-3630QM CPU 2.4 GHz. 


quantity 

PCI 

PC2 

PC3 

PC4 

Qi 

18.14 

19.51 

20.64 

35.39 

n num 
" i 

13.53 

21.40 

26.78 

35.56 


Table 6 : Salt content in any small brick in the treated case (PC). We reported the salt content 
qi and the numerical values q r / urn expressed in mg/cm 3 . 


Let us define the average porosity in i-th brick as: 


1 


—t / n(z, t)dz for i = 1, 2, 3,4, (51) 

K Ja t 

then for completeness, we report in the next Table [7] the average porosity obtained numerically 
n”“ m , i = 1,..., 4 for the four bricks obtained discretizing (51) with the trapezoidal quadrature 
rule, both in the not treated (NT) and treated (PC) case. 
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brick 1 


brick 2 



brick 3 




brick 4 



Figure 9: Experiment 2. Comparison between the profile of c s without and in presence of 
PC-10“ 6 M. 


porosity 

^num 
n l 

^num 
ri 2 

num 

n 3 

num 

'*4 

NT 

0.2355 

0.2324 

0.2280 

0.2251 

PC 

0.2331 

0.2343 

0.2264 

0.2233 


Table 7: Porosity in any small brick. We reported the average value for the porosity n™ um ,i = 
1,..., 4 for the not treated and the treated bricks. 

6 Conclusions 

We developed a mathematical model to describe the action of crystallization inhibitors into a 
porous stone. This simple model is able to capture the main features of the inhibitor from 
experiments carried out on a set of commercially available bricks. According to the current 
knowledge, the model describes the action of inhibitors through two coefficients: crystalliza¬ 
tion rate, I\ s , taking into account nucleation, and the specific volume 7 , taking into account 
the crystal habit modification. From the calibration of the mathematical model described in 
section [3j we found out that the action of phosphocitrate (PC) increases the crystallization 
rate and decreases the crystal specific volume. This means that, although crystals form faster 
in the presence of the inhibitor, nevertheless they occupy a smaller volume, thus lowering the 
development of tensile stresses, and, on the other hand, ensuring the hydraulic continuity into 
the porous stones. In the future, we will repeat the same study varying the materials and with 
more detailed experiments in order to test and improve our mathematical model. Our aim is 
to end up with a sound simulation tool to investigate crystallization modifier. 
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